Evolution of pairing from weak to strong coupling on a honeycomb lattice 
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We study the evolution of the pairing from weak to strong coupling on a honeycomb lattice by 
Quantum Monte Carlo. We show numerical evidence of the BCS-BEC crossover as the coupling 
strength increases on a honeycomb lattice with small fermi surface by measuring a wide range 
of observables: double occupancy, spin susceptibility, local pair correlation, and kinetic energy. 
Although at low energy, the model sustains Dirac fermions, we do not find significant qualitative 
difference in the BCS-BEC crossover as compared to those with an extended Fermi surface, except 
at weak coupling, BCS regime. 
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I. INTRODUCTION 

It has long been known that the pairing formed from 
an attractive coupling has a smooth crossover between 
the weak coupling and the strong coupling^^. In the 
weak coupling limit, singlet pairs are formed around 
the fermi surface, according to the BCS theory. In the 
strong coupling limit, local bound pairs can be formed, 
and these "preformed pairs" condense as the tempera- 
ture is further lowered where the Bose-Einstein condcn- 
sation(BEC) occurs. The interest on this crossover has 
been revitalized^^ii^^, mainly due to the quest of 
understanding the pseudogap phase in the high temper- 
ature superconductors. 

Recently, condensed matter systems sustain on 
fermions with linear dispersion, typical examples are hon- 
eycomb lattice models and nodal fermions for d-wave 
superconductors, have generated huge surge of inten- 
sive studies. These models possess substantial differ- 
ences from models with extended Fermi surface such 
as models on square lattice. In particular, it has been 
suggested that the quantum phase transition (QPT) be- 
tween the metallic phase and the degenerate charge den- 
sity wave/pairing phase at half- filling in the attractive 
Hubbard model (AHM) on honeycomb lattice is related 
to its BCS-BEC crossover away from half-fillingJ 1 . This 
certainly does not happen on the square lattice, in which 
the flat Fermi surface at half-filling renders the Umklapp 
scattering becoming the dominant channel, its BCS-BEC 
crossover is not related to any QPT through tuning the 
attractive coupling—. In the honeycomb lattice, the den- 
sity of state is zero at half-filling, therefore any instability 
from the band structure is weakened, and strong cou- 
pling is needed to induce ordering. It can be shown that 
all the short range interactions are irrelevant. In order 
to tackle the strong coupling problem, besides breaking 
the symmetry by mean field ansatz, we choose Quantum 
Monte Carlo method in this work to study the BCS-BEC 
crossover in the honeycomb lattice. 

Various studies£ ' 13 i 14 ' 15 have been devoted to the BCS- 



BEC crossover of the AHM on a square lattice. The ob- 
jective of this work is to study how do the linear disper- 
sion, and the aforementioned QPT at half-filling affect 
the BCS-BEC crossover of the slightly doped system. 

Our main finding can be summarized as follow. At 
the weak coupling, BCS-like regime, pseudogap phenom- 
ena are observed, however we expect that it is mainly 
due to the band structure of honeycomb lattice, rather 
than the bound pair formation. At the intermediate cou- 
pling, crossover regime, we can identify two temperature 
scales, the high temperature one where the performed 
pair formed with associated pseudogap phenomena; and 
the low temperature one where the system enters the 
pairing phase. At strong coupling, BEC-like regime, the 
electrons form pairs at high temperature and condense 
as hard core bosons at low temperature. However, we 
do not find distinctive feature compares to the square 
lattice, except at the weak coupling regime where the 
band structure dominates the quasi-particle dispersion. 
Further interpretations of the QMC results are next pre- 
sented by applying the mean field (MF) approximation 
to lattice models and continuum model for fermions with 
linear dispersion. 

II. MODEL AND METHOD 

The AHM in honeycomb lattice reads 

H=-t c ta c 3° - U^riifnn -/x^rw, (1) 

<i, i 7>,tr i ia 

where Ci a (cf a ) annihilates (creates) a particle with spin 
a at site i, denotes the nearest-neighbor lattice 

sites i and j, t is the hopping matrix element, U is the 
on-site attractive interaction, and [i is the chemical po- 
tential. In the following we set t = 1 as the energy 
scale of the system, all the observable are in units of 
t. The bare electronic (U = limit) dispersion is given 

by e k = ±Js + 2cos(V3k y ) + 4cos(V3k y /2)cos(3k x /2), 
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and the band width W is 6. At half filling this is linear 
around the Fermi points. Keeping only the low energy 
excitations, in the first quantized form the wave function 
follows the 2D Weyl equation for masslcss chiral Dirac 
fermions, vpa ■ V , 5(r) = Ei£(r), where a — (a x ,a- y ) are 
the Pauli matrices and vf — 3/2 is the Fermi velocity. 
This description in term of Dirac fermions is not exact 
away from half filling. Nevertheless, the linear dispersion 
can be a good approximation below the van Hove singu- 
larities at filling n = 1 ± 1/4. For this reason, we choose 
n = 0.88 for our calculations using determinant quantum 
Monte Carlo(DQMC)±&. 

The DQM C 16 i 17 is a Hamiltonian based approach. 
The Hamiltonian H in the partition function Z = 
Tr exp(— ftH ) is expressed in the real space via the Trot- 
ter decomposition and Hubbard-Stratonovich (HS) trans- 
formation. The only systematic error is from discretizing 
the imaginary-time ft into M slices of At = ft/M in the 
Trotter decomposition. The HS transformation replaces 
the on-site interactions in the attractive Hubbard model 
by HS fields coupled to the charge. The summation over 
the HS fields is treated by Monte Carlo procedure. The 
calculations are proceeded on a N — 72 sites honeycomb 
lattice, the actual lattice for the simulation is shown in 
the Fig. [T] Since the attractive Hubbard coupling does 
not have minus-sign problem, a wide range of tempera- 
tures and couplings can be studied. 
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FIG. 1: Sketch of a 72 sites honeycomb lattice. The red and 
green solid circles are the lattice points in the honeycomb 
lattice. The red solid circles also represent the underlying 
triangular lattice. T\ and T2 are the real space translational 
vectors. 



III. QMC RESULTS 

One of the clear signals indicating the formation of 
bound pairs at strong coupling is the formation of spin 
gap. At weak coupling, we expect fermion quasi-particle 
character to remain at high temperature, for which the 



spin susceptibility increases as the temperature is low- 
ered. On the other hand, the strong coupling limit is 
manifested by the decrease of the spin susceptibility as 
the temperature is lowered, due to the formation of the 
gap which leads to the reduction in the spectral function 
at low frequency. 

We first show the spin susceptibility x(q, at fre- 
quency oj — 0, and momentum q = (0, 0) in Fig. [2j 
where we also show the spin susceptibility from RPA cal- 
culation for comparison. %(0, 0) is suppressed for all cou- 
plings, as can be inferred simply from the RPA formu- 
lation, where xrpa(0,0) = Xo(0,0)/(1 + U X o(0,0)). At 
weak coupling x(0, 0) increases as the temperature is low- 
ered as expected for a fermion quasi-particle description, 
however it bends downward before it goes upward again 
as the temperature is lowered further. This two peak 
structure of x(0, 0) associated with the formation of the 
pscudogap has been found in the dynamical mean field 
theory study^. However, in the honeycomb lattice, the 
apparent pseudogap phenomena indicated by this struc- 
ture of x(0, 0) already exist in the weak coupling regime, 
below the strong coupling regime where the "preformed 
pair" phenomena occur. Therefore, we believe that it is 
derived from the particular dispersion relation of honey- 
comb lattice, where the density of state is small around 
the doped Fermi surface. 

On the other hand, in the strong coupling regime, 
x(0, 0) vanishes quickly as the bound pairs are formed 
and spin gap equals the binding energy needed to break 
the pair. In the weak coupling regime, the QMC re- 
sults behave similarly as compared to the RPA results. 
When the interaction is increased to around W/2, the 
QMC results evolve in the opposite direction as com- 
pared to the RPA results and drop sharply at low tem- 
perature, whereas the RPA results at low temperature 
limit do not change qualitatively when U increases. This 
signals that the system enters the phase in which elec- 
trons form bound pairs, and the spin excitations start 
to be gapped^. The pairing phase cannot be reached 
by summing the ladder diagrams within the RPA. For 
strong coupling (U « W), the suppression of x(0,0) be- 
comes smooth and appears at high temperature. This 
effect reflects the fact that the bound pairs are already 
formed at high temperature^. The temperature where 
deviations appear between the QMC results and the RPA 
results is an indication of the formation of local singlet 
pair, which can be interpreted as the energy scale where 
the fermion quasi-particle description is not valid for any 
lower temperature. 

We then probe the pairing directly by considering 
the pairing correlation function for local pairing, Pq = 

■m S;,i( c !+/,T c !+/a c ^l Ci .T + h - c -)- Tne onl y instability 
is pairing, in this incommensurate doped case (rules out 
CDW order). We expect Pq to increase as temperature 
is lowered for all coupling strengths. One of the most 
representative characteristics of local pairs is that they 
are distributed uniformly in space and condense around 
zero momentum as bosons when temperature is lowered. 
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FIG. 2: Uniform spin susceptibility x(0, 0) (left), and pairing 
correlation function Po — Po (right) as a function of temper- 
ature for a range of interaction strength at n = 0.88. 



This is manifested by the rapid increase of local pair cor- 
relation as shown in Fig. [2j note that the single particle 
contribution Po has been subtracted from Po to empha- 
size the vertex contribution of pairing 18 . The conden- 
sation of the bosonic local pairs for the strong coupling 
case (U = 5, 6) can be observed from the rapid increase 
of P — Po with the decrease of temperature. In contrast, 
the pairs formed around the quasi-particle Fermi surface 
in the weak coupling regime only bring a slight increase 
in P - P . 
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FIG. 3: The kinetic energy (left), and double occupancy (right) 
as a function of temperature for a range of interaction strength 
at n = 0.88. 

We show the kinetic energy, Ek = 

{-tIN)Y, H) ,M,° c ^) in Fi s- 13 In the weak 

coupling regime, its temperature dependence is similar 
to the free fermion case. When we increase the inter- 
action to the crossover regime (U « 3 — 4), qualitative 
change already happened in the high temperature, 
where the gain in the kinetic energy is much slower than 



the free fermion case. Moving into the strong coupling 
regime, fermions begin to form bound pairs at high 
temperature and only lose little kinetic energy. When 
temperature further decreases, the local pairs in the 
system condense and hence Ek drops sharply; 14 .. 

A good indicator to measure the local pair formation in 
the BEC state is the double occupancy (n^ni), see Fig. 
[3] We find that {n^n\) increases as the temperatures 
decrease. However, it reaches a local maximum at cer- 
tain temperature. This can be understood as the change 
of the kinetic energy which destabilizes the double occu- 
pancy. This behavior of (n^n{) are in accord with the 
fact that the local maximum coincides with the temper- 
ature where the kinetic energy drops most sharply. At 
very low temperature, the bosonic on-site pairs begin to 
dominate, (njnj.) increases again and should saturate at 
n/2 for strong couplings. 

After elaborating the evidence of BCS-BEC crossover, 
we put those observables from DQMC together and iden- 
tify the temperature scales for different U . In Fig. [4] we 
show the results represented for weak (U = 1), interme- 
diate (U — 3,4), and strong (U — 6) couplings. 

At U — 1, x(0:0) QMC result does not deviate from 
the RPA result. The local pair correlation does not de- 
velop, and (n^n^) is small even at low temperature, which 
shows that the pairing correlation is weak. The critical 
temperature, T c , for the Kosterlitz-Thouless transition 
into the pairing phase is below the lowest temperatures 
we studied. 

At U = 3, x(0,0) local pair correlations begin to in- 
crease at T c 0.2. At almost the same temperature 
the QMC result begins to deviate from the RPA result. 
These imply the developments in both spin and pairing 
correlations. The system shows BCS-likc pairing effect 
from the instability of the fcrmi surface. However, there 
is no true phase coherence at any finite temperature as 
that in the BCS theory. Nevertheless, at this coupling 
strength, the pairing is still rather weak, due to the small 
density of state around the Fermi energy. 

At U = 4 the system displays two temperature scales. 
The first one is T* at high temperature around T « 0.8, 
this could be associated with the pseudogap phase. At 
this temperature, x(0, 0) from QMC result reaches its 
maximum and begins to deviate from the RPA result. 
In addition (n^jii) also reaches the first plateau at high 
temperature. These signal that electrons bound pairs 
start to develop, spin gap is formed and the quasi-particle 
description is broken below this temperature. We es- 
timate the critical temperature for the condensation of 
bound pairs, T c ss 0.3. Below this temperature, the local 
pair correlation Po — Po grows quickly and x(0, 0) drops 
sharply; (n^ni) reaches its low temperature maximum 
and saturates. 

At U = 6, the system is at the strong coupling limit, 
where U reaches the band width W, there is only one 
temperature scale in the system, T c ss 0.5, within the 
temperature range we studied. x(fii 0) reaches its maxi- 
mum at very high temperature and decreases smoothly, 
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FIG. 4: Double occupancy (n^n^), uniform spin susceptibility 
x(0, 0), and pairing correlation Po -Po as a function of tem- 
perature for different coupling strength at (n) — 0.88 filling. 
The magenta shadow regions are used to mark the energy 
scale. 

which suggests that pair formation begins at a very high 
temperature, above the temperature range we studied. 
Below T c , Po — -Po increase quickly, and (n^ni) tends 
to n/2 at zero temperature. These suggest that the 
bound pairs undergo a Kosterlitz-Thouless transition, 
which manifests a BEC-like scenario. 

From the above numerically exact DQMC data, we 
show clearly that there is a qualitative change from weak 
to strong coupling at finite temperature. This should 
correspond to the true BCS-BEC crossover at zero tem- 
perature. However, we find that the results for the hon- 
eycomb lattice have no drastic qualitative difference as 
compared to that of the square lattice^. Certainly, the 
band structure alters the quantitative values of the cou- 



pling for the crossover. However, the BCS-BEC crossover 
on a doped honeycomb lattice models exists at U « 3 — 4 
where the linear dispersion approximation for the free 
fermions is not valid. 

IV. DISCUSSIONS AND CONCLUSIONS 

With the progress of the techniques of optical lat- 
tices and the fabrication of single layer graphene, the 
BCS-BEC crossover on a honeycomb lattice and Dirac 
fermions is not only an important problem itself, but also 
has broad experimental and theoretical interests with 
other topics under intensive studies. 

The atom-atom interaction in ultracold fcrmionic 
atoms in a optical trap can be tuned by magnetic field 
Feshbach resonance. The honeycomb lattice can pos- 
sibly be realized by optical trapi^. This may provide 
a direct way to study experimentally the BCS-BEC 
crossover problem with linear dispersion.— . In addi- 
tion, the superconducting phase of graphene via the at- 
traction from phonons and plasmons has been discussed 
recentl y 20 ' 21 i 22 ' 23 i 24 . Although it is unlikely to generate 
strong attraction from phonon coupling in graphene, our 
results suggest that even at weak coupling regime, non- 
trivial temperature dependence of spin susceptibility may 
occur in the superconducting phase from local Holstcin 
phonon coupling. 

In conclusion, we have presented extensive results from 
DQMC which confirm the BCS-BEC crossover for the 
doped (n=0.88) AHM on a honeycomb lattice. In con- 
trast to the systems with extended fermi surface, there is 
an enhancement of pseudogap property revealed from the 
double peak structure in the spin susceptibility at weak 
coupling due to the peculiar density of state of honey- 
comb lattice. Apart from this, the BCS-BEC crossover 
does not show prominent difference between square lat- 
tice and honeycomb lattice for the parameters and system 
size we study. 
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